# TABLES A9-A12

# Author: Kasia Nalewajko
# First created: 6 February 2023
# Replicated: 8 June 2024

rm(list = ls())

# LOAD PACKAGES -----------------------------------------------------------

if (!require("dplyr")) install.packages("dplyr")
if (!require("tidyr")) install.packages("tidyr")
if (!require("stringr")) install.packages("stringr")
if (!require("lmtest")) install.packages("lmtest")
if (!require("sandwich")) install.packages("sandwich")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("stargazer")) install.packages("stargazer")

# LOAD DATA ---------------------------------------------------------------

load("./00 SUBMITTED/00 APSR final/04 replication_files/01 data/RANfrance.Rda")

# RUN MODELS ---------------------------------------------------------------

hiding_bi <- lm(data = frran_hh,
                formula = hiding ~ type)
hiding_deptFE <- lm(data = frran_hh,
                    formula = hiding ~ type + as.factor(dept_code))
hiding_districtFE <- lm(data = frran_hh,
                        formula = hiding ~ type + as.factor(deppct))

organising_evasion_bi <- lm(data = frran_hh,
                            formula = organising_evasion ~ type)
organising_evasion_deptFE <- lm(data = frran_hh,
                                formula = organising_evasion ~ type + as.factor(dept_code))
organising_evasion_districtFE <- lm(data = frran_hh,
                                    formula = organising_evasion ~ type + as.factor(deppct))

finding_shelter_bi <- lm(data = frran_hh,
                         formula = finding_shelter ~ type)
finding_shelter_deptFE <- lm(data = frran_hh,
                             formula = finding_shelter ~ type + as.factor(dept_code))
finding_shelter_districtFE <- lm(data = frran_hh,
                                 formula = finding_shelter ~ type + as.factor(deppct))

forging_documents_bi <- lm(data = frran_hh,
                           formula = forging_documents ~ type)
forging_documents_deptFE <- lm(data = frran_hh,
                               formula = forging_documents ~ type + as.factor(dept_code))
forging_documents_districtFE <- lm(data = frran_hh,
                                   formula = forging_documents ~ type + as.factor(deppct))

# CALCULATE ROBUST STANDARD ERRORS AND EXPORT THE MODELS

rob.fit1 <- coeftest(hiding_bi, function(x) vcovHC(x, type="HC0"))
rob.fit2 <- coeftest(hiding_deptFE, function(x) vcovHC(x, type="HC0"))
rob.fit3 <- coeftest(hiding_districtFE, function(x) vcovHC(x, type="HC0"))

cov1 <- sandwich::vcovHC(hiding_bi, type = "HC1")
cov2 <- sandwich::vcovHC(hiding_deptFE, type = "HC1")
cov3 <- sandwich::vcovHC(hiding_districtFE, type = "HC1")

robust_se1 <- sqrt(diag(cov1))
robust_se2 <- sqrt(diag(cov2))
robust_se3 <- sqrt(diag(cov3))

stargazer::stargazer(list(hiding_bi, hiding_deptFE, hiding_districtFE),
                     se = list(robust_se1, robust_se2, robust_se3),
                     title="Linear Probability Models of likelihood of performing specific rescue activities, by rescuer profile",
                     align=TRUE,
                     model.numbers=TRUE,
                     keep = c("type"),
                     dep.var.labels="Dependent variable: Hiding (binary)",
                     type = "latex")
                     
rob.fit4 <- coeftest(organising_evasion_bi, function(x) vcovHC(x, type="HC0"))
rob.fit5 <- coeftest(organising_evasion_deptFE, function(x) vcovHC(x, type="HC0"))
rob.fit6 <- coeftest(organising_evasion_districtFE, function(x) vcovHC(x, type="HC0"))

cov4 <- sandwich::vcovHC(organising_evasion_bi, type = "HC1")
cov5 <- sandwich::vcovHC(organising_evasion_deptFE, type = "HC1")
cov6 <- sandwich::vcovHC(organising_evasion_districtFE, type = "HC1")

robust_se4 <- sqrt(diag(cov4))
robust_se5 <- sqrt(diag(cov5))
robust_se6 <- sqrt(diag(cov6))

stargazer::stargazer(list(organising_evasion_bi, organising_evasion_deptFE, organising_evasion_districtFE), # ols1, ols2,
                     se = list(robust_se4, robust_se5, robust_se6),
                     title="Linear Probability Models of likelihood of performing specific rescue activities, by rescuer profile",
                     align=TRUE,
                     model.numbers=TRUE,
                     keep = c("type"),
                     dep.var.labels="Dependent variable: Organising Evasion (binary)",
                     type = "latex")

rob.fit7 <- coeftest(finding_shelter_bi, function(x) vcovHC(x, type="HC0"))
rob.fit8 <- coeftest(finding_shelter_deptFE, function(x) vcovHC(x, type="HC0"))
rob.fit9 <- coeftest(finding_shelter_districtFE, function(x) vcovHC(x, type="HC0"))

cov7 <- sandwich::vcovHC(finding_shelter_bi, type = "HC1")
cov8 <- sandwich::vcovHC(finding_shelter_deptFE, type = "HC1")
cov9 <- sandwich::vcovHC(finding_shelter_districtFE, type = "HC1")

robust_se7 <- sqrt(diag(cov7))
robust_se8 <- sqrt(diag(cov8))
robust_se9 <- sqrt(diag(cov9))

stargazer::stargazer(list(finding_shelter_bi, finding_shelter_deptFE, finding_shelter_districtFE), # ols1, ols2,
                     se = list(robust_se7, robust_se8, robust_se9),
                     title="Linear Probability Models of likelihood of performing specific rescue activities, by rescuer profile",
                     align=TRUE,
                     model.numbers=TRUE,
                     keep = c("type"),
                     dep.var.labels="Dependent variable: Finding shelter (binary)",
                     type = "latex")

rob.fit10 <- coeftest(forging_documents_bi, function(x) vcovHC(x, type="HC0"))
rob.fit11 <- coeftest(forging_documents_deptFE, function(x) vcovHC(x, type="HC0"))
rob.fit12 <- coeftest(forging_documents_districtFE, function(x) vcovHC(x, type="HC0"))

cov10 <- sandwich::vcovHC(forging_documents_bi, type = "HC1")
cov11 <- sandwich::vcovHC(forging_documents_deptFE, type = "HC1")
cov12 <- sandwich::vcovHC(forging_documents_districtFE, type = "HC1")

robust_se10 <- sqrt(diag(cov10))
robust_se11 <- sqrt(diag(cov11))
robust_se12 <- sqrt(diag(cov12))

stargazer::stargazer(list(forging_documents_bi, forging_documents_deptFE, forging_documents_districtFE), # ols1, ols2,
                     se = list(robust_se7, robust_se8, robust_se9),
                     title="Linear Probability Models of likelihood of performing specific rescue activities, by rescuer profile",
                     align=TRUE,
                     model.numbers=TRUE,
                     keep = c("type"),
                     dep.var.labels="Dependent variable: Forging Documents (binary)",
                     type = "latex")
                     
